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Space-based gravitational-wave (GW) detectors, such as LISA or a similar ESA-led mission, will 
offer unique opportunities to test general relativity. We study the bounds that space-based detectors 
could realistically place on the graviton Compton wavelength At, = h/{rngc) by observing multiple 
inspiralling black hole (BH) binaries. We show that while observations of individual inspirals will 
yield mean bounds Ag ~ 3 x 10^* km, the combined bound from observing ~ 50 events in a two-year 



mission is about ten times better: Ag ~ 3 x 10 km (r 



4 X 10" 



eV). The bound improves 



faster than the square root of the number of observed events, because typically a few sources provide 
constraints as much as three times better than the mean. This result is only mildly dependent on 
details of BH formation and detector characteristics. The bound achievable in practice should 
be one order of magnitude better than this figure (and hence almost competitive with the static, 
model-dependent bounds from gravitational effects on cosmological scales), because our calculations 
ignore the merger /ringdown portion of the waveform. The observation that an ensemble of events 
can sensibly improve the bounds that individual binaries set on Ag applies to any theory whose 
deviations from general relativity are parametrized by a set of global parameters. 



The formulation of gravitational theories with nonzero 
mass for the graviton that are consistent with cosmologi- 
cal observations is an important open problem. Attempts 
to construct such theories led to well-known concep- 
tual difficulties, such as the so-called van Dam-Vcltman- 
Zakharov (vDV-Z) discontinuity [i|-y|, due to the fact 
that the helicity-0 component of the graviton does not 
decouple from matter when the putative mass of the 



Our analysis will show that, in hierarchical models of 
massive BH formation [y, 0] , the bound on Ag from GW 
observations of a population of merger events is about 
one order of magnitude better than the mean bound from 
GW observations of individual mergers. 



graviton mg 



0. To circumvent pathologies related 



to the vDV-Z discontinuity, various versions of Lorentz- 
violating massive graviton theories have been proposed 
in recent years [^. Massive graviton signatures in the 
CMB and possible constraints on rrig from cosmological 
observations are an active area of research (see e.g. Q)- 
In this paper we are interested in hypothetical massive 
graviton theories as "straw men" for alternative theo- 
ries of gravity in which the propagation speed of grav- 
ity differs from that of electromagnetic waves, leading to 
a modified dispersion relation. Therefore we will adopt 
a phenomenological point of view and ask the following 
question: i/the graviton mass were nonzero, what upper 
bounds on rrig could we set by gravitational-wave (GW) 
observations of inspiralling compact binaries with future 
space-based detectors? Using Ag = h/{mgc), upper lim- 
its on the graviton mass rrig (in cV) can be expressed as 
lower limits on its Compton wavelength Ag (in km): 



Current bounds 


Ag[fcm] 


mg[eV] 


Reference 


Binary pulsars* 


1.6 X 10^° 


7.6 X 10-20 


[8j 


Solar System^ 


2.8 X 10^2 


4.4 X 10-22 


\g,M 


Clusters^ 


6.2 X 10"/io 


2.0 X lO-28/i-i 


[1^ 


Weak lensingt 


1.8 X 10^2 


6.9 X lO-^'^ 


[12] 


Proposed bounds 


Ag[fcm] 


m,[eV] 


Reference 


Pulsar timing* 


4.1 X 10^^ 


3.0 X 10-23 


[13] 


White dwarfs* 


1.4 X 10" 


8.8 X 10-24 


[14] 


EM counterparts* 


^Ql5 _ -^Ql6 


10-24 - 10-25 


[15] 



Ag[km] X mg[cV] = 1.24 x 10' 



(1) 
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TABLE I. Graviton mass bounds. For proposed methods we 
quote the best achievable bounds. In the notation of the main 
text, a dagger ^ denotes static bounds; a number sign *, dy- 
namical bounds; an asterisk * , bounds that could be achieved 
by comparing GW and electromagnetic observations. 



To put our results in context, in Table |T] we present a 
non- exhaustive summary of current and proposed bounds 
on Ag that do not rely solely on GW observations. 

These bounds can be roughly divided into three classes. 
i) Static bounds. If the graviton has nonzero mass, 
it is reasonable to expect that the Newtonian gravita- 
tional potential will be modified to the Yukawa form 



in the nonradiative near zone of any body of mass M: 
V{r) = {GM/r) exp {-r/Xg) . Talmadge et al. @ inves- 
tigated deviations from Kepler's third law for the inner 
planets of the Solar System. By translating the uncer- 
tainties in these measurements in terms of accelerations 
of test bodies resulting from a Yukawa potential, Will 
|10| found that the strongest bound on Ag comes from the 
very nearly Kcplerian orbit of Mars; Xg > 2.8 x 10^'^ km. 

Some bounds on mg quoted by the Particle Data 
Group (PDG) [ifll similarly assume Yukawa corrections 
to the Newtonian potential in the weak-field limit. The 
strongest bounds are naturally obtained from observa- 
tions on the largest (cosmological) scales. As early as 
1974, Goldhaber and Nieto assumed that the graviton 
mass would produce a Yukawa-type correction to the 
standard Newtonian potential and argued that the ev- 
idence for bound clusters and tidal interactions between 
galaxies should imply a range for gravity at least as large 
as a few Mpc, so that Xg > 6.2 x lO^^ft-o km [11|. How- 
ever, the need to include dark matter to explain galactic 
rotation curves indicates that there are complications in 
the nature of gravity on those scales that are not neces- 
sarily well characterized by a Yukawa-type potential. An 
even stronger bound A^ > 1.8 x 10^^ km comes from weak 
gravitational Icnsing, because no distortions are observed 
in the measured values of the variance of the power spec- 
trum [l^. Because of uncertainties in the amount and 
dynamics of dark matter in the Universe (and in the ab- 
sence of a consistent massive graviton theory compati- 
ble with cosmology) these bounds should be regarded as 
mo del- dependent and viewed with some caution. 
ii) Dynamical bounds. All bounds listed in part i) are 
static, in the sense that they do not probe features related 
to the propagation of the gravitational interaction when 
rUg ^ 0. The best available dynamical bounds come from 
the (indirect) observations of GWs from binary pulsars, 
that are in excellent agreement with general relativity 
pTf . Finn and Sutton [^ observed that the consistency of 
the orbital decay rates of binary pulsars PSR B1534+12 
and B1913+16 with general relativistic predictions yields 
a "dynamical" bound Ag > 1.6 x 10^° km. A recent idea 
is to set bounds on rrig using Pulsar Timing Arrays [13| . 
The best bound achievable in the near future would be 
Ag ~ 4.1 X 10"'^'^ km, but this figure could worsen by an 
order of magnitude depending on the number of pulsars 
used for the test, timing accuracy and observation time. 
iii) Comparisons of gravitational and electromag- 
netic observations. If nig ^ 0, the modified dis- 
persion relation for GWs would result in different ar- 
rival times of GWs and electromagnetic waves emitted 
by the same astrophysical source. Cutler ei al. [1J| 
proposed to correlate electromagnetic observations and 
future space-based GW observations of white dwarf bi- 
naries. The best bound that could be obtained in this 
way is Ag ^ 1.4 x 10"'^"' km, but realistic bounds would 
probably be worse by about one order of magnitude. 

Kocsis et al. proposed to correlate LISA observations 
of GWs from massive BH binary mergers with their pos- 



sible electromagnetic counterparts [l5J . An intrinsic limi- 
tation of this method is related to timing uncertainties in 
the GW burst, which are comparable to the dynamical 
timescale for the coalescing binary during merger, and 
can be estimated as the inverse of the orbital frequency 
at the innermost stable circular orbit. For binaries in 
the mass range Al = 10^ — 10^ at redshift z = 1, this 
uncertainty leads to a best bound Ag ~ 10^^ — 10^® km, 
worse than bounds coming from GW observations alone 
(as we will see below). Systematic, model-dependent un- 
certainties in the electromagnetic counterpart will further 
weaken graviton mass bounds achievable in this way. 



I. BOUNDS ON \g FROM INDIVIDUAL 
GRAVITATIONAL- WAVE OBSERVATIONS 

Tight dynamical bounds on the graviton mass can be 
obtained using GW observations in space. This is due 
to two reasons: (1) the larger mass of observable BH bi- 
naries, as compared to ground-based GW observations; 
(2) the statistical increase in the bound that would re- 
sult from observing populations of individually resolved 
binaries. To illustrate point (1), in this section we review 
existing work on massive graviton bounds from GW ob- 
servations of individual binaries with LISA. In section |TT] 
we will present the first attempt to quantify the statistical 
improvement of bounds on Ag that could be achievable 
in reality by observing several events in the context of 
hierarchical models of massive BH formation. 

Will [13] first pointed out that interesting bounds on 
Ag could come from a careful monitoring of the phase 
of GWs emitted by binaries of compact objects, such as 
BHs and/or neutron stars. In hypothetical massive gravi- 
ton theories, the GW damping formulae and dispersion 
relation would be modified. As a consequence, the GW 
phasing ^gr(/) would acquire an additional term: 



*mg(/)-*gr(/)-/3sM^/)" 



(2) 



where /3g = 7t'^DAI/[{1+z)X'1] and D is a, distance param- 
eter, similar to (but not quite the same as) the luminosity 
distance Dl (here and below we assume a ACDM model 
with Ho = 70 km-s-^Mpc-i, ^m = 0.3, Qm = 0.7). 

Using the restricted post-Newtonian approximation 
and neglecting binary spins. Will estimated that stellar- 
mass binary inspirals to be observed with LIGO would 
yield a bound Ag ~ 5 x 10^^ km, only slightly better than 
Solar System bounds. For an "optimal" LISA system, i.e. 
an equal-mass BH binary of total mass M = 2 x 10^ Mq 
at Dl = 3 Gpc (z ~ 0.5), the bound would be four orders 
of magnitude better: Ag ~ 5 x 10^^ km ~ 1.6 kpc. 

These initial estimates were refined in various papers. 
Will and Yunes [l8| showed that bounds from binary BH 
inspirals at Dl = 3 Gpc would range between 10^^ km 
and 5 x lO^^ km for M in the range 10^ - 10^ Mq. 
They found that the bound on Ag is proportional to vL 
(where L is the LISA armlength) and to (LISA acceler- 
ation noise)"'^/^, and that it scales in the following way: 



XgOC 



D \'/' Xii/12 



il + z)DL 



1/4 .1/3 



"^0 ^0 



(3) 



where M = (77117712)^/^ A/^^^^ is the "chhp mass", Sq (m 
Hz~^) sets the scale of the noise power spectral density 
(PSD) of the detector, and /o is a characteristic "knee" 
frequency where the PSD has a minimum. Distance de- 
pendence is weak because the effect of the massive gravi- 
ton and measurement errors both grow with distance. 

Berti et al. [13 performed a more detailed Monte 
Carlo calculation of LISA's parameter estimation capa- 
bilities. They still considered restricted post-Newtonian 
inspiral waveforms, but they used Cutler's model [201 to 
take into account the motion of the detector (all previ- 
ous analyses assumed a static LISA constellation). In 
this paper we will do the same. In Cutler's model, 
LISA can be seen as one (two) independent "LIGO- 
likc" Michclson interferometers depending on whether 4 
(5/6, respectively) laser links are available between the 
three satellites forming the constellation. For the "op- 
timal binary" (a nonspinning, equal-mass binary with 
M = 2 X 10*^ Mq at Di = 3 Gpc), Monte Carlo averages 
over the binary position and orientation yield a bound 
Xg = 5.0 X 10^'^(3.7 X 10^*^) km when using two (one) 
Michelsons, in excellent agreement with Will's original 
results. For masses in the range M = 2 x 10''— 2 x lO"^ Mq , 
the mean bound using two Michelsons at Dl = 3 Gpc is 
in the range ~ 5 x 10^^ - 7 x lO^^ km. 

Later work used mostly Cutler's model and investi- 
gated the effect of more complex inspiral waveforms. 
Arun and Will 2l[ studied the effect of higher harmonics 
and PN amplitude corrections. They found that bounds 
on Xg improve by a factor of a few for systems of total 
mass > 10^ Mq, and that this improvement is more sig- 
nificant for binaries with large mass ratios. If one ignores 
spin precession, adding spins to the waveforms generally 
introduces degeneracies between the binary parameters, 
degrading parameter estimation accuracy [l9[. Stavridis 
and Will [22[ showed that modulations induced by spin 
precession can break these degeneracies. For the "opti- 
mal binary" , for example, the average bound including 
spin precession is Ag ^ 5 x 10^^ km, basically the same as 
in the case of nonspinning binaries. Yagi and Tanaka |23| 
included both spin precession and eccentricity. Perform- 
ing simulations for a (10^ -|- 10^) Mq binary at 3 Gpc, 
they found an average bound Xg = 3.1 x 10^^ km. This 
is again consistent with [1^ to within a factor 2. 

More recently, Kcppel and Ajith [2^ revisited this 
problem using phenomenological waveforms for nonspin- 
ning BH binaries that include the merger/ringdown 
phase. Their analysis of LISA bounds is similar to the 
original work by Will [10[ in that it ignores the motion of 
the detector, using an effective non-sky-averaged noise 
PSD. This simplification should not affect their main 
conclusions: (i) for the "optimal binary" , the bound 
on Xg improves by about one order of magnitude, up 



to ^ 4 X 10^^ km; (ii) a comparable bound is obtained 
also for binaries of larger mass, up to M ^ 10^ Af©. 
At Dl = 3 Gpc, the best bound using the full merger 
is Xg ^ 5.9 X 10^^ km for M = 4.8 x lO"^ Mq, to be 
compared with a best bound Xg = 6.3 x 10^^ km for 
M = 1.9 X 10^ Mq if we consider only inspiral waves. A 
similar order-of-magnitude improvement is expected for 
Earth-based detectors. 

Del Pozzo et al. [2^ recently revisited graviton mass 
bounds using Bayesian inference (see also [26l - [29f for sim- 
ilar work including other possible alternatives to general 
relativity). They focused on ground-based observations 
of sources within 150 Mpc and found results consistent 
with Will's original analysis [i3|: for example, their Fig. 6 
shows that binary observations with Advanced LIGO 
would yield typical bounds ^ few x 10^^ km. 

In this study we consider quasicircular, nonspinning, 
restricted post-Newtonian inspiral waveforms. We con- 
sider an observation time of two years. We take into ac- 
count weak lensing errors on the redshift following [30|, 
and we compute the individual graviton mass bounds 
by generalizing the Fisher matrix formalism described in 
[3l| , which takes into account correlations with the other 
waveform parameters. The results of Refs. [22, [2J] sug- 
gest that our bounds will be (i) very close to the bounds 
we would obtain for spinning, precessing inspirals, and 
(ii) about one order of magnitude worse than the bounds 
achievable if we used "full" merger waveforms. 



II. BOUNDS ON \g FROM 

GRAVITATIONAL- WAVE OBSERVATIONS OF 

MASSIVE BLACK HOLE POPULATIONS 



All studies of massive graviton bounds so far analyzed 
isolated systems, such as Will's "optimal binary" of mass 
M ~ 2 X 10^ Mq at z ^ 0.5. Unfortunately, hierarchical 
models of massive BH formation and evolution predict 
that typical systems observable by space-based interfer- 
ometers would have masses smaller than this, and be 
located at redshift z ~ 4 or higher (see e.g. [y, 0| and 
Fig. 2 of |31|). The main uncertainties in these models 
concern the seeding mechanism and the role of accretion 
in BH growth [3l|, |33 • In this work we ignore spins in the 
gravitational waveforms. Accretion mostly influences the 
spin magnitude J33| , so we will focus on the role of seed- 
ing. Following work by the LISA Parameter Estimation 
Taskforce [3J], we will consider two "extreme" scenarios: 
small seeds, efficient accretion (SE) and large seeds, ef- 
ficient accretion (LE). These models being extreme, we 
expect that our results should bracket the constraints 
that would be obtained using other population models. 

Here we consider the "Classic LISA" design along with 
two different designs for the proposed ESA-led space- 
based detector, that we will call by the working name 
of "New LISA" . "Classic LISA" consists of three space- 
craft forming an equilateral triangle with laser power 
P = 2 W, telescope diameter d = 0.4 m and armlength 
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FIG. 1. Non sky-averaged noise power spectral density 
S^'^^U) for New LISA C2 (black, solid line), New LISA C5 
(blue, dash-dotted line) and Classic LISA (red, dashed line). 




FIG. 2. Distribution of bounds on the graviton Compton 
wavelength for individual observations. We consider 1000 re- 
alizations of the massive BH population and three different 
detector designs: New LISA C2 (black, thick line). New LISA 
C5 (blue, medium line) and Classic LISA (red, dashed line). 



Mean (median) of individual events (10^ km) 



Detector 



Classic LISA 
New LISA C2 
New LISA C5 



SE, 1 Mich. LE, 1 Mich. SE, 2 Mich. LE, 2 Mich 



4.26(2.60) 6.83(5.77) 
3.03(2.44) 3.62(3.27) 
3.41(2.53) 4.63(4.13) 



Mean (median) of combined bound (10 km) 



Detector 



Classic LISA 
New LISA C2 
New LISA C5 



SE, 1 Mich. LE, 1 Mich. 



4.93(4.87) 5.67(5.59) 
2.29(2.25) 2.73(2.71) 
3.10(3.07) 3.64(3.62) 



L = 5 X 10® m, trailing 20° behind the Earth at an in- 
clination of 60° with respect to the ecliptic. The authors 
are members of a Science Performance Task Force that 
is considering several different LISA-likc configurations 
with different characteristics and sensitivities. The con- 
figurations that we call New LISA C2 (C5, respectively) 
consist of three spacecraft forming an equilateral triangle 
with armlength L = 10® m (L = 2 x 10® m), laser power 
P = 2 W and telescope diameter d = 0.4 m (d = 0.28 m). 
"New LISA" should be deployed 10° behind the Earth, 
gradually drifting to ~ 25° behind the Earth in 5 years. 

The non sky- averaged noise power spectral densities 
for all three configurations are shown in Fig. [TJ they 
are related to the sky-averaged power spectral density by 
S^^^if) = ^S^^{f) (see [li] for a discussion of sky av- 
eraging). These curves include galactic confusion noise, 
estimated using methods similar to [3J| (which in turn 
was based on [35|). In our study we consider the noise 
power spectral density to be infinite below a cutoff fre- 
quency /cutoff = 10~^ Hz. As shown in [l9|, bounds on 
Xg drop significantly at masses > 2 x 10^ Mq if the noise massive BH population, in units of 10^*^ km. 
cannot be trusted below lO^"' Hz. This assumption has 
a mild effect on our results, because most binaries in our 
models have mass lower than this. 

For each model we consider 1000 realizations of the 
Universe. Each of these realizations typically produces 
~ 30 — 50 events observable with signal-to-noise ratio 
(SNR) larger than 8. The distribution of bounds on Xg 
resulting from individual observations is shown in Fig. [51 

The top half of Tabic |lT] shows that the mean bound 
over individual observations is ^ 3 x 10^^ km for New 
LISA C2, and only slightly better for the other designs. 
This conclusion is quite robust, in the sense that numbers 
vary only mildly for different seeding mechanisms and 
different detector characteristics. 



4.87(2.72) 9.13(7.72) 
3.60(2.80) 4.76(4.29) 
4.02(2.84) 6.15(5.48) 



SE, 2 Mich. LE, 2 Mich. 



6.51(6.45) 7.50(7.37) 
3.09(3.04) 3.66(3.64) 
4.16(4.12) 4.85(4.82) 



TABLE II. Top: mean (in parentheses: median) bound on \g 
for different BH formation models, using one or two detectors, 
in units of 10^^ km. Bottom: mean (in parentheses: median) 
of the combined bound on A^ over 1000 realizations of the 



In most alternative theories, deviations from general 



relativity can be parametrized by one or more global pa- 
rameters (such as Xg) which are the same for every sys- 
tem. It is natural to expect that one can obtain better 
constraints on these parameters, as well as other univer- 
sal constants, by combining multiple observations (see 
e.g. [2l,[2l[2i,[3i-l3i|). Assuming that estimates for in- 
dividual sources are independent and Gaussian posteriors 
for each source, consistent with the Fisher matrix approx- 
imation, the width a^ of the combined posterior on l/Xg 
is given by a^^ ~ Si '^r^' where ct| is the width of the 
posterior for the ith source. The bound on Xg can thus be 
obtained by adding the individual bounds in quadrature. 
The results are shown in Fig. |3] and in the bottom half of 
Table [ITl The combined bound obtained from the whole 
BH population is about one order of magnitude better 
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FIG. 3. Distribution of combined bounds over 1000 realiza- 
tions of the MBH population. Linestyles are as in Fig. [2] 



than the average bound obtained from typical observa- 
tions. A rough estimate would suggest that N identical 
sources should provide a bound ^ VN times better than 
the bound from a single source. Our combined bound 
is typically about 3 times better than the bound from 
the best event, but the median bound is typically an or- 
der of magnitude viforse than the best, and hence ^ 30 
times worse than the combined bound. A typical realisa- 
tion has ~ 50 events, so our analysis shows that we can 
beat the yN extrapolation from the median bound by 
a considerable margin. If 5/6 links (two Michelsons) are 
available instead of 4 links (one Michelson), the bound 
typically improves by a factor ~ v^- 



III. CONCLUSIONS AND OUTLOOK. 

We assessed the capability of future space-based in- 
terferometers, such as "Classic LISA" and the proposed 
ESA-led "New LISA" , to constrain the mass of the gravi- 
ton by combining observations of a population of massive 
BH binaries. We found that: (1) by using a population 
of merging BH binaries we can obtain a bound on Xg that 
is ^ 10 times better than the mean bound on individual 
observations; (2) quite independently of the detector's 
design and of details of the massive BH formation mod- 
els, the combined bound from inspiral observations will 
be Ag ~ 3 X 10^^ km. This figure is likely to under- 
estimate the bound achievable in practice by about one 
order of magnitude, as we have ignored the merger and 
ringdown portion of the waveform [2J] , but further work 
is required to confirm this expectation. 

In conclusion, space-based observations of a popula- 
tion of merging BHs should set bounds in the range 
Xg E [2 X 10^^ , 10^^ km] on the graviton Compton wave- 
length, depending on details of the detector and on the 
specific waveform model used to set the bounds. This is 
comparable to the (static and model-dependent) bounds 
from cosmological-scale observations quoted in Table H] 
but it is very different in nature, because gravitational 
radiation tests the dynamical regime of Einstein's gen- 
eral relativity. 
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